#include <iostream>
#include <fstream>
#include <cmath>
#include "/home/mjarvis/misc/nr/ran0.cpp"
#include "/home/mjarvis/misc/nr/zbrent.cpp"
#include "Int.h"
#include "dbg.h"
using std::endl;

std::ostream* dbgout=new std::ofstream("makewwa.debug");
bool XDEBUG = false;

const char* cmbfile = "CMBQuintchainNew.dat";

double glob_omega_m;
double glob_w0;
double glob_wa;

double cmb_dist_integrand(double a)
{
  double lna = log(a);
  double h = glob_omega_m * exp(-3.*lna) + 
    (1.-glob_omega_m) * exp( -3.*(1.+glob_w0+glob_wa)*lna - 3.*glob_wa*(1.-a) );
  h = sqrt(h);
  double asq = a*a;
  return 1./(h*asq);
}

double cmb_dist(double omega_m, double w0, double wa)
// d_a = a_cmb int( 1/H(a) a^(-2) da , a=a_cmb..1)
// H(a) = H0 sqrt( omega_m a^(-3) +
//                 omega_DE exp( 3 int(1+w(a') dln(a'),a'=a..1) ) )
// With w(a) = w0 + wa(1-a), this becomes:
// H(a) = H0 sqrt( omega_m a^(-3) +
//                 omega_DE exp( -3(1+w0+wa)ln(a) - 3wa(1-a) ) )
// Since we are looking for changes in this quantity for different
// parameters, I actually calculate d_a*H0/a_cmb
{
  const double acmb = 1./1090.;

  glob_omega_m = omega_m;
  glob_w0 = w0;
  glob_wa = wa;
  integ::IntRegion<double> reg(acmb,1.);
  double dist = integ::int1d(std::ptr_fun(cmb_dist_integrand),reg);
  //dbg<<"dist("<<omega_m<<","<<w0<<","<<wa<<") = "<<dist<<endl;
  return dist;
}
  
const double w_min = -10.;
const double w_max = 4.;
const double wa_min = -4.;
const double wa_max = 8.;

double glob_da0;

double daroot(double w)
{
  // Want to solve for da(w) = da0
  // so root equation is da(w) - da0
  return cmb_dist(glob_omega_m,w,glob_wa)-glob_da0;
}

double FindW(double omega_m,double wa,double da)
{
  glob_da0 = da;
  glob_omega_m = omega_m;
  glob_wa = wa;
  double neww = zbrent(daroot,w_min,w_max,1.e-4);
  dbg<<"da("<<omega_m<<","<<neww<<","<<wa<<") = ";
  dbg<<cmb_dist(omega_m,neww,wa)<<endl;
  return neww;
}

double range_rand(double min, double max)
{
  return min + (max-min)*ran0();
}

int main()
{
    std::ifstream fin(cmbfile);
    if (!fin) myerror("opening cmb file: ",cmbfile);

    double ombh2,omch2,theta,n,Z,w0,As,h,ll;
    while (fin >> ombh2 >> omch2 >> theta >> n >> Z >> w0 >> As >> h >> ll) {
        double omega_m = (omch2+ombh2)/(h*h);
        double da = cmb_dist(omega_m,w0,0.0);
        dbg<<"da("<<omega_m<<","<<w0<<",0.0) = "<<da<<endl;
        double wa = range_rand(wa_min,wa_max);
        double neww = FindW(omega_m,wa,da);
#if 0
        double newwplussig = FindW(omega_m,wa,da-0.1);
        double newwminussig = FindW(omega_m,wa,da+0.1);
        double sig_w = std::abs(newwplussig - newwminussig)/2.;
        dbg<<"wa = "<<wa<<"  sig_w = "<<sig_w<<endl;
        // Use sigma_w for the likelihood to account for the phase space
        // volume at that value of wa.
        double like = sig_w;
#endif
        double like =  1.;

        std::cout << ombh2 <<"  "<< omch2 <<"  "<<theta<<"  "<<n<<"  ";
        std::cout << wa<<"  "<<neww<<"  "<<As<<"  "<<h<<"  "<<like<<std::endl;
    }
}
